The effect of a dominant shrub on annuals plants in a desert system changes along a gradient of precipitation.
Alex Filazzola
May 2017
Ephedra regional gradient (ERG) analyses
Abstract
An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.
Hypothesis
We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.
Methods
load GAMs Gradient
Year<-c("2016","2016","2016","2016","2016","2016","2016","2017","2017","2017","2017","2017","2017","2017")
Site<-c("Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch","Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch")
gamsq <-c(7992.36,6789.76,7903.21,5715.36,7779.24,7482.25,7447.69,6872.41,5882.89,7157.16,5595.04,7293.16,7106.49,6955.56)
gams.richard <-c(89.4,82.4,88.9,75.6,88.2,86.5,86.3,82.9,76.7,84.6,74.8,85.4,84.3,83.4)
gams <- c(89.6,71,84.7,62.7,88.3,80.6,80.7,63.6,51.4,73.4,49.5,69.4,68.3,65.8)
gams.data <- data.frame(Year,Site,gamsq,gams)Weather for growing season
Climate patterns within study
season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)
### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days
season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)## Site frost.days
## 1 Barstow 24.725275
## 2 Cuyama 29.670330
## 3 HeartofMojave 25.824176
## 4 PanocheHills 10.439560
## 5 SheepholeValley 6.043956
## 6 Tecopa 23.626374
## 7 TejonRanch 50.549451season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)## Site frost.days
## 1 Barstow 17.679558
## 2 Cuyama 19.889503
## 3 HeartofMojave 16.574586
## 4 PanocheHills 14.917127
## 5 SheepholeValley 1.785714
## 6 Tecopa 25.966851
## 7 TejonRanch 35.911602## Both years had comparable number of frost days
## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}
season1.frost <- season1 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 10
## 2 Cuyama 10
## 3 HeartofMojave 12
## 4 PanocheHills 5
## 5 SheepholeValley 7
## 6 Tecopa 11
## 7 TejonRanch 14season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 7
## 4 PanocheHills 6
## 5 SheepholeValley 3
## 6 Tecopa 7
## 7 TejonRanch 13## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)## Site frost.days avg.min.temp
## 1 Barstow 12.396694 5.750413
## 2 Cuyama 11.570248 3.352893
## 3 HeartofMojave 16.528926 4.542149
## 4 PanocheHills 2.479339 6.777686
## 5 SheepholeValley 2.479339 9.394167
## 6 Tecopa 11.570248 6.342149
## 7 TejonRanch 33.884298 1.438017season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)## Site frost.days avg.min.temp
## 1 Barstow 11.666667 6.052500
## 2 Cuyama 16.666667 3.624167
## 3 HeartofMojave 8.333333 5.637838
## 4 PanocheHills 8.333333 6.065833
## 5 SheepholeValley 0.000000 9.386916
## 6 Tecopa 20.833333 5.938333
## 7 TejonRanch 28.333333 2.535833season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 4
## 4 PanocheHills 2
## 5 SheepholeValley 2
## 6 Tecopa 4
## 7 TejonRanch 10season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 3
## 4 PanocheHills 3
## 5 SheepholeValley 2
## 6 Tecopa 5
## 7 TejonRanch 10Differences from Nutrient content
# nutrient data for each site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
## join aridity with nutrients
nutrients.climate <- merge(nutrients,subset(gams.data, Year==2016), by=c("Site"))
## Nitrogen difference between shrub and sites
m.nit <- aov(log(N) ~ Site * microsite, data=nutrients)
summary(m.nit)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 63.04 10.51 13.241 2.86e-09 ***
## microsite 1 40.16 40.16 50.614 2.24e-09 ***
## Site:microsite 6 4.96 0.83 1.042 0.408
## Residuals 56 44.43 0.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.nit, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(N) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 1.514873 1.088317 1.941429 0## Potassium difference between shrub and sites
m.pot <- aov(log(K) ~ Site * microsite, data=nutrients)
summary(m.pot)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 15.897 2.650 22.013 4.09e-13 ***
## microsite 1 4.445 4.445 36.934 1.14e-07 ***
## Site:microsite 6 1.605 0.268 2.223 0.0541 .
## Residuals 56 6.740 0.120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.pot, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(K) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 0.5040081 0.3378755 0.6701406 1e-07## Phosphorus difference between shrub and sites
m.pho <- aov(log(P) ~ Site * microsite, data=nutrients)
summary(m.pho)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 23.163 3.861 14.33 8.10e-10 ***
## microsite 1 10.546 10.546 39.15 5.79e-08 ***
## Site:microsite 6 3.394 0.566 2.10 0.0676 .
## Residuals 56 15.085 0.269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.pho, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(P) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 0.7762734 0.5277317 1.024815 1e-07nutrient.mean <- nutrients.climate %>% group_by(Site, gams) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)
ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=nitrogen), color="#E69F00", size=3) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#E69F00") + geom_jitter(aes(x=gams, y=potassium/20), color="#56B4E9", size=3) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=potassium/20), color="#56B4E9") + geom_jitter(aes(x=gams, y=phosphorus), color="#009E73", size=3) +ylab("soil nutrient content")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("rainfall continentality")+
annotate("segment", x = 72, xend = 74, y = 25, yend = 25, colour = "#56B4E9", size=2) + ## add custom legend
annotate("text", x = 65, y = 25, label = "Nitrogen", size=5, hjust=0)+
annotate("segment", x = 72, xend = 74, y = 23, yend = 23, colour = "#E69F00", size=2) + ## add custom legend
annotate("text", x = 65, y = 23, label = "Potassium", size=5, hjust=0)+
annotate("segment", x = 72, xend = 74, y = 21, yend = 21, colour = "#009E73", size=2) + ## add custom legend
annotate("text", x = 65, y = 21, label = "Phosphorus", size=5, hjust=0)
m1 <- lm(nitrogen~poly(gams,2), data=nutrient.mean)
summary(m1)##
## Call:
## lm(formula = nitrogen ~ poly(gams, 2), data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## 1.4326 3.1541 -1.6134 -1.2346 0.2699 -0.2336 -1.7750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4326 0.8325 6.526 0.00285 **
## poly(gams, 2)1 10.6460 2.2026 4.833 0.00844 **
## poly(gams, 2)2 9.0253 2.2026 4.098 0.01488 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.203 on 4 degrees of freedom
## Multiple R-squared: 0.9094, Adjusted R-squared: 0.8641
## F-statistic: 20.08 on 2 and 4 DF, p-value: 0.008208m2 <- lm(phosphorus~gams, data=nutrient.mean)
summary(m2)##
## Call:
## lm(formula = phosphorus ~ gams, data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## 3.9609 1.3929 -1.9844 0.5018 -0.6999 -5.4853 2.3140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.6031 11.4769 2.057 0.0948 .
## gams -0.1929 0.1432 -1.347 0.2357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.399 on 5 degrees of freedom
## Multiple R-squared: 0.2664, Adjusted R-squared: 0.1196
## F-statistic: 1.815 on 1 and 5 DF, p-value: 0.2357m3 <- lm(potassium~poly(gams,2), data=nutrient.mean)
summary(m3)##
## Call:
## lm(formula = potassium ~ poly(gams, 2), data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## 8.0453 -13.5433 0.9712 4.7180 -15.7101 25.3691 -9.8502
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 186.557 6.705 27.824 9.93e-06 ***
## poly(gams, 2)1 -244.894 17.740 -13.805 0.000160 ***
## poly(gams, 2)2 159.200 17.740 8.974 0.000853 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.74 on 4 degrees of freedom
## Multiple R-squared: 0.9855, Adjusted R-squared: 0.9782
## F-statistic: 135.6 on 2 and 4 DF, p-value: 0.0002114Native non-native comparisons
### annual plant community repsonse
spp.data <- read.csv("Data/ERG.communitydata.csv")
spp.data[is.na(spp.data)] <- 0
## load community data
community <- read.csv("Data/ERG.communitydata.csv")
community[is.na(community)] <- 0
## load species list
spp.list <- read.csv("Data/ERG.specieslist.csv")
## convert data to long format
status <- gather(community, species, abundance, 13:53)
names(spp.list)[1] <- "species" ## rename species column for merge
## combine native-non-native status with data set
status <- merge(status, spp.list, by="species")
## summarize per plot native and non-natives, convert back to wide format
mean.status <- status %>% group_by(Year, Site, Microsite, Rep, status) %>% summarize(abundance=sum(abundance))
mean.status <- spread(mean.status, status, abundance)
status.data <- mean.status
## calculate site means for native and non-native
mean.status <- mean.status %>% group_by(Year, Site, Microsite) %>% summarize(native=mean(native), non.native=mean(non.native))Phylogenetics
library(picante)
library(ape)
library(brranching)
library(taxize)
## load species list to create tree
spp.list <- read.csv("Data/ERG.specieslist.csv")
## create a combined species name column and drop erinoginum because not to species level
taxa <- paste(spp.list$Genus, spp.list$Species.name)
taxa <- taxa[taxa!="Erinoginum spp"]
## create tree of phylogeny
tree <- phylomatic(taxa=taxa, get = 'POST')
## calculate branch lengths using Grafen's method
## Grafen, A. (1989) The phylogenetic regression. Philosophical Transactions of the Royal society of London. Series B. Biological Sciences, 326, 119–157.
tree2 <- compute.brlen(tree, method="Grafen")
## view tree and outputs to verify branches
plot(tree2)
tree2##
## Phylogenetic tree with 39 tips and 29 internal nodes.
##
## Tip labels:
## cryptantha_barbigera, cryptantha_intermedia, phacelia_crenulata, phacelia_tanacetifolia, pectocarya_spp, amsinckia_grandiflora, ...
## Node labels:
## poales_to_asterales, eudicots, , , asterids, ericales_to_asterales, ...
##
## Rooted; includes branch lengths.## format community to only have site and microsite
comm <- community %>% group_by(Year, Site, Microsite) %>% summarise_if(is.numeric, funs(sum(., na.rm=T)))
comm <- comm[,c(1:3,13:53)] ## extract species abundances only
comm <- data.frame(comm) ## convert to dataframe
comm[,"site.micro"] <- paste(comm$Site,comm$Microsite,as.character(comm$Year)) ## create site by microsite column
rownames(comm) <- comm[,length(comm)] ## replace row names with site by microsite
comm[,c("Year","Site","Microsite","Erinoginum.spp","Boraginaceae.sp.","site.micro")] <- NULL ## drop all columns but ones for analysis
## match species name format
spp.list[,"spp.name"] <- paste(spp.list$Genus, spp.list$Species.name)
colnames(comm) <- spp.list$spp.name[match(colnames(comm), spp.list$Species.shorthand)]
## match species format
names(comm) <- gsub(" ", "_", names(comm), fixed=T) ## underscores instead of spaces
names(comm) <- gsub("__", "_", names(comm), fixed=T) ## replace double underscores with one
names(comm) <- tolower(names(comm)) ## lower case names
## mean phylogenetic distance
mpd.data <- mpd(comm, cophenetic(tree2), abundance.weighted=T)
## format dataframe
mpd.data <- data.frame(Year=c(rep(2016,14),rep(2017,14)),Microsite=c("open","shrub"),site=rownames(comm), mpd=mpd.data)
mpd.data[,"Site"] <- gsub(" open", "", mpd.data[,"site"])
mpd.data[,"Site"] <- gsub(" shrub", "", mpd.data[,"Site"])
mpd.data[,"Site"] <- gsub(" 2016", "", mpd.data[,"Site"])
mpd.data[,"Site"] <- gsub(" 2017", "", mpd.data[,"Site"])
mpd.data[,"site"] <- NULL
mpd.data[is.na(mpd.data)] <- 0 ## barstow no plants
## getspecies richness for sites
spp.rich <- community %>% group_by(Year, Site, Microsite) %>% summarize(richness=mean(Richness))
spp.rich <- data.frame(spp.rich)
## join data
mpd.rich <- merge(mpd.data, spp.rich, by=c("Site","Microsite","Year"))Community analysis
### indicator species analysislibrary(indicspecies)
## indicator species analysis
indval <- multipatt(community[,13:53], community$Microsite, control=how(nperm=999))
summary(indval)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 41
## Selected number of species: 15
## Number of species associated to 1 group: 15
##
## List of species associated to each combination:
##
## Group open #sps. 8
## stat p.value
## E.cicutarium 0.558 0.001 ***
## L.gracilis 0.427 0.001 ***
## A.wrangelianus 0.201 0.001 ***
## C.brevipes 0.171 0.036 *
## A.lentiginosus 0.149 0.027 *
## C.rigida 0.149 0.015 *
## A.grandiflora 0.129 0.021 *
## Boraginaceae.sp. 0.120 0.027 *
##
## Group shrub #sps. 7
## stat p.value
## C.fremontii 0.379 0.001 ***
## M.glabrata 0.277 0.001 ***
## P.tanacetifolia 0.267 0.001 ***
## A.tessellata 0.199 0.003 **
## B.nigra 0.179 0.017 *
## E.nitidum 0.172 0.004 **
## M.affinis 0.129 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1site.names <- as.character(unique(community$Site))
community.site <- subset(community, Site==site.names[7])
## species accumilation curves
accum1 <- specaccum(subset(community.site, Microsite=="shrub", select=13:53), method="random", permutations=1000)
accum2 <- specaccum(subset(community.site, Microsite=="open", select=13:53), method="random", permutations=1000)
par(mar=c(4.5,4.5,0.5,0.5))
plot(accum2, col="#E69F00",ci.col="#E69F0090", xlab="Sampling effort", cex.axis=1.5, cex.lab=1.8, ylab="Species richness", lwd=2, ylim=c(0,16))
plot(accum1,add=T , col="#56B4E9", ci.col="#56B4E990", lwd=2)
micros <- community$Microsite[1:120]
indval <- multipatt(community.site[,13:53], micros, control=how(nperm=999))
summary(indval)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 41
## Selected number of species: 4
## Number of species associated to 1 group: 4
##
## List of species associated to each combination:
##
## Group open #sps. 1
## stat p.value
## Boraginaceae.sp. 0.316 0.031 *
##
## Group shrub #sps. 3
## stat p.value
## A.tessellata 0.531 0.001 ***
## B.nigra 0.475 0.011 *
## M.affinis 0.343 0.042 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cluster analysis
library("cluster")
library("factoextra")
library("magrittr")
comm.sum <- community %>% group_by(Microsite, Site) %>% summarize_if(is.numeric, funs(sum))
comm.sum <- data.frame(comm.sum)
rownames(comm.sum) <- paste(comm.sum$Microsite,comm.sum$Site)
## coreelation matrix
### Distance measures
res.dist <- get_dist(comm.sum[,13:53], stand = TRUE, method = "euclidean")
fviz_dist(res.dist, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
## dendrogram
## transform data
comm.trans <- decostand(comm.sum[,13:53], "hell")
## clean data for ordination
## see distribution of spp
boxplot(comm.trans, xaxt="n")
labs <- colnames(comm.trans)
text(cex=0.8, x=1:41-1, y=-0.12, labs, xpd=TRUE, srt=45)
## remove spp with only one instancve
comm.trans <- comm.trans[,!colSums(comm.trans)==apply(comm.trans, 2, max)]
## replace outliers with mean
avg.max <- function(x) {
y = max(x)
avg = mean(x)
x[x==y] <- avg
return(x)
}
#
# ## dot chart to find outliers (remove or Avg)
# #dotchart(comm.trans$E.cicutarium) ## Erodium
# comm.trans[,"E.cicutarium"] <- avg.max(comm.trans$E.cicutarium) ## removed outlier
# #dotchart(comm.trans$A.wrangelianus) ## Chilean trefoil
# comm.trans[,"A.wrangelianus"] <- NULL ## removed entire species
# #dotchart(comm.trans$A.lentiginosus) ##
# comm.trans[,"A.lentiginosus"] <- avg.max(comm.trans$A.lentiginosus) ## removed outlier
# #dotchart(comm.trans$H.vulgare) ## Hordeum
# comm.trans[,"H.vulgare"] <- NULL ## removed entire species
# #dotchart(comm.trans$Pectocarya.spp) ##
# comm.trans[,"Pectocarya.spp"] <- avg.max(comm.trans$Pectocarya.spp) ## removed outlier
# #dotchart(comm.trans$L.arizonicus) ## Lupin
# comm.trans[,"L.arizonicus"] <- NULL ## removed entire species
# #dotchart(comm.trans$M.bellioides) ##
# comm.trans[,"M.bellioides"] <- NULL ## removed entire species
# #dotchart(comm.trans$B.nigra) ## Mustard
# comm.trans[,"B.nigra"] <- NULL ## removed entire species
# #dotchart(comm.trans$N.demissum) ## Sand Matt
# comm.trans[,"N.demissum"] <- avg.max(comm.trans$N.demissum) ## removed outlier
# #dotchart(comm.trans$L.gracilis) ##
# comm.trans[,"L.gracilis"] <- avg.max(comm.trans$L.gracilis) ## removed outlier
# #dotchart(comm.trans$A.grandiflora) ## Amsinckia
# comm.trans[,"A.grandiflora"] <- NULL ## removed entire species
# #dotchart(comm.trans$B.diandrus) ## Bromus diandrus
# comm.trans[,"B.diandrus"] <- NULL ## removed entire species
# #dotchart(comm.trans$M.affinis) ##
# comm.trans[,"M.affinis"] <- NULL ## removed entire species
# #dotchart(comm.trans$P.crenulata) ##
# comm.trans[,"P.crenulata"] <- NULL ## removed entire species
# #dotchart(comm.trans$Erinoginum.spp) # buckwheat
# comm.trans[,"Erinoginum.spp"] <- avg.max(comm.trans$Erinoginum.spp) ## removed outlier
# #dotchart(comm.trans$C.lasiophyllus) #
# comm.trans[,"C.lasiophyllus"] <- NULL ## removed entire species
#
#
# library(usdm)
library(corrplot)
#
# ## Check for collinearity
# vifcor(comm.trans)
#
# ## drop Claviformis because of collinear problems in P.insularis
# comm.trans[,"C.claviformis"] <- NULL ## removed entire species
# comm.trans[,"L.decorus"] <- NULL ## removed entire species
# comm.trans[,"E.glyptosperma"] <- NULL ## removed entire species
# comm.trans[,"E.wallacei"] <- NULL ## removed entire species
# ## outliers removed
# cor.dat <- cor(comm.trans[,1:20])
# corrplot(cor.dat, method="number")
## outliers left
cor.dat <- cor(comm.trans[,1:35])
corrplot(cor.dat, method="number")
## calculate distances
dis <- vegdist(comm.trans, method="bray")
## cluster analysis
clus <- hclust(dis, "ward.D2")
## cluster colours
colz <- c("#3e4444", "#86af49", "#405d27","#c1946a")
fviz_dend(clus, k = 4, # Cut in four groups
cex = 0.7, cex.axis=2, # label size
k_colors = colz,
color_labels_by_k = TRUE, # color labels by groups
rect = TRUE # Add rectangle around groups
) 
grp <- cutree(clus, 4)
## Interpretation of classes
## get soil moisture from site/microsite
## Soil moisture
smc <- read.csv("Data/ERG.phytometer.census.csv")
smc.end <- subset(smc, Census=="end")
smc.avg <- smc.end %>% group_by(Site, Microsite) %>% summarize(mean.smc=mean(swc, na.rm=T))
smc.avg <- data.frame(smc.avg)
smc.avg[,"grp"] <- as.factor(c(3,3,2,1,4,4,1,1,4,4,3,3,2,2))
boxplot(mean.smc ~ grp, data=smc.avg)
## CA or PCA
dca1 <- decorana(comm.trans) ## length of gradient >2 & determine relative differences in community composition
## conduct ordination
ord <- cca(comm.trans)
summary(ord)##
## Call:
## cca(X = comm.trans)
##
## Partitioning of mean squared contingency coefficient:
## Inertia Proportion
## Total 1.957 1
## Unconstrained 1.957 1
##
## Eigenvalues, and their contribution to the mean squared contingency coefficient
##
## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.6995 0.3254 0.2379 0.1838 0.15873 0.13238 0.08478
## Proportion Explained 0.3574 0.1663 0.1216 0.0939 0.08112 0.06765 0.04332
## Cumulative Proportion 0.3574 0.5238 0.6453 0.7392 0.82033 0.88799 0.93131
## CA8 CA9 CA10 CA11 CA12 CA13
## Eigenvalue 0.03921 0.03664 0.02865 0.01793 0.009198 0.002791
## Proportion Explained 0.02004 0.01872 0.01464 0.00916 0.004700 0.001430
## Cumulative Proportion 0.95135 0.97007 0.98471 0.99387 0.998570 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CA1 CA2 CA3 CA4 CA5 CA6
## B.rubens 1.5386 -0.58468 -0.16290 -0.300662 0.032522 0.025724
## E.cicutarium 0.8449 0.03763 0.23090 0.198152 -0.159488 0.651322
## A.wrangelianus 1.3890 -0.25748 0.30331 0.348959 -0.598941 2.555549
## S.barbatus -0.1070 0.76275 -0.18049 -0.290922 -0.130899 -0.038848
## C.barbigera -0.4246 0.26354 -0.19786 -0.001936 -0.207406 -0.040565
## C.claviformis -0.7806 -0.30349 0.26357 -0.159678 0.017537 -0.066644
## A.lentiginosus -0.5098 -0.20775 1.06481 -0.537354 -0.027420 0.228868
## L.arizonicus -0.9562 -0.72017 -1.36251 0.652822 0.109890 0.084549
## C.rigida -0.4416 0.54538 0.42769 -0.658482 -0.279382 -0.006550
## H.vulgare 1.8619 -1.42186 -0.31281 -0.499631 -1.199167 -1.279525
## P.tanacetifolia 1.2704 -1.21579 0.01562 -0.057382 -1.407379 -2.124678
## Pectocarya.spp -0.5912 -0.11588 0.63679 -0.119010 0.001115 -0.082389
## D.capitatum 0.9839 0.46497 0.62533 2.049871 -1.030658 -0.663977
## C.brevipes -0.4353 0.06706 0.96613 0.262212 -0.152629 -0.026021
## L.decorus -0.6477 -0.46992 1.37440 0.057320 0.110230 -0.051132
## P.insularis -0.8224 -0.52038 0.32355 -0.025486 0.107794 0.004392
## M.glabrata -0.7516 -0.39205 -1.07168 0.394740 0.162632 0.126873
## C.fremontii -0.8106 -0.42600 -0.37872 0.185389 0.048897 0.100343
## C.intermedia -0.3248 0.09038 0.73133 -0.398634 0.373265 -0.212349
## M.bellioides -0.7993 -0.67433 1.37262 -0.249245 0.199870 -0.099398
## B.nigra 1.0327 0.63789 -0.18882 -0.087370 2.109966 -0.503093
## E.wallacei -0.8720 -0.70776 0.35312 0.085296 0.206398 0.050102
## N.demissum -0.7986 -0.49732 0.95631 -0.260716 0.113486 -0.047266
## L.gracilis 0.8954 0.74787 0.46052 1.545242 -0.046187 -0.344203
## A.tessellata 1.0882 0.53061 -0.14406 0.038462 2.053034 -0.346446
## A.grandiflora 1.3468 -0.17135 0.34349 0.499486 -0.606150 2.413659
## P.secunda 1.5501 -0.68523 -0.16117 -0.158957 -0.236351 -0.496685
## E.glyptosperma -0.9120 -0.79133 -0.34064 0.388338 0.224171 0.069830
## B.diandrus 1.0706 0.57659 -0.21693 -0.121689 2.267130 -0.528020
## M.affinis 1.0631 0.58866 -0.21140 -0.114931 2.236184 -0.523112
## L.sparsiflorus -0.7748 -0.72143 -1.25994 0.568275 -0.045061 -0.114310
## P.crenulata -0.9387 -0.90518 -0.59025 0.567872 0.304319 0.120311
## Erinoginum.spp -0.7619 0.05689 -0.88747 -0.053460 -0.240620 -0.141904
## E.nitidum -0.8606 -0.66499 0.04005 0.178949 0.154017 -0.029840
## C.lasiophyllus -0.4208 1.12726 -0.69866 -0.670429 -0.652207 0.110872
##
##
## Site scores (weighted averages of species scores)
##
## CA1 CA2 CA3 CA4 CA5 CA6
## open Barstow -0.2084 1.4664 -0.48332 -0.9412 -0.85976 0.35431
## open Cuyama 0.6663 1.2156 0.99043 2.9233 -0.72223 0.12889
## open HeartofMojave -0.8113 -0.6255 1.67608 -0.4331 0.24940 0.07784
## open PanocheHills 1.5677 -0.6215 0.13350 -0.2873 -0.56847 3.15526
## open SheepholeValley -0.9135 -0.5324 -1.20730 0.4543 -0.06212 -0.06157
## open Tecopa -0.4560 1.2568 0.06852 -1.1593 -0.65414 -0.34438
## open TejonRanch 0.7591 1.0798 0.01388 0.1600 0.97701 -0.32340
## shrub Barstow -0.4288 1.0054 -0.81161 -0.4759 -0.66361 0.26214
## shrub Cuyama 1.7446 -1.3327 -0.24905 -0.0420 -1.76930 -2.56278
## shrub HeartofMojave -0.7912 -0.7075 1.16636 -0.1243 0.16621 -0.21986
## shrub PanocheHills 2.0376 -1.5554 -0.40828 -1.1849 -0.34542 0.64209
## shrub SheepholeValley -1.0243 -1.0200 -1.61034 0.9698 0.38452 0.31785
## shrub Tecopa -0.3991 1.4576 -0.39266 -1.1976 -0.62130 -0.29896
## shrub TejonRanch 1.1930 0.3788 -0.30764 -0.2324 2.77412 -0.60843## calculate priority
spp.priority <- colSums(comm.trans)
## plot ordination
par(mar=c(4.5,4.5,0.5,0.5))
plot(ord, type="n")
ordiellipse(ord, grp, lty = 2, col = "grey80", draw="polygon", alpha=150)
orditorp(ord, display = "species", cex = 0.7, col = "darkorange3", priority=spp.priority, air=0.8)
orditorp(ord, display = "sites", cex = 0.7, col = "darkslateblue", air=0.1)
##collect environmental variables for site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
nutrients.mean <- nutrients %>% group_by(microsite, Site) %>% summarise_all(funs(mean))
nutrients.vars <- data.frame(nutrients.mean)
## summarize daily means across sites
HOBO <- read.csv("Data/ERG.logger.data.csv")
HOBO.data <- HOBO %>% group_by(Site, Microsite, Year, Month, Day) %>% summarize(Temp=mean(Temp),RH=mean(RH))
means <- HOBO.data %>% group_by(Microsite, Site) %>% summarize(Temp.=mean(Temp, na.rm=T),rh=mean(RH, na.rm=T),temp.se=se(Temp),rh.se=se(RH))
means <- data.frame(means)
## soil moisture
smc.early <- subset(smc, Census=="emergence")
smc.avg <- smc.early %>% group_by(Microsite, Site) %>% summarize(SMC=mean(swc, na.rm=T))
envs <- data.frame(swc=smc.avg[,"SMC"],nutrients.vars[,c("N","P","K")], means[,c("Temp.","rh")])
## Check for collinearity
cor(envs)## SMC N P K Temp. rh
## SMC 1.0000000 -0.31692385 0.44854712 0.81694831 -0.6372354 0.9185204
## N -0.3169238 1.00000000 0.05879919 -0.05275422 0.4167329 -0.4367433
## P 0.4485471 0.05879919 1.00000000 0.64885814 -0.5303896 0.5342666
## K 0.8169483 -0.05275422 0.64885814 1.00000000 -0.5769100 0.7750256
## Temp. -0.6372354 0.41673289 -0.53038959 -0.57690998 1.0000000 -0.8650006
## rh 0.9185204 -0.43674328 0.53426662 0.77502564 -0.8650006 1.0000000envs[,"K"] <- NULL ## drop potassium for being correlated
envs[,"rh"] <- NULL ## drop humidity for being correlated
ord.env <- envfit(ord, envs)
plot(ord.env)
Short-term vs long-term effects
Repeated of analyses separating long- and short- term effects
Ordination including aridity and microsite (pCCA)
Partial Constrained Correspondence Analysis for community data
Compare Ambient Community - mixed model results
community.arid <- merge(community, gams.data, by=c("Year","Site"))
mean.spp <- community.arid %>% group_by(Microsite, Site, gams, Year) %>% summarize(abd=mean(Abundance),rich=mean(Richness), bio=mean(Biomass))
mean.status2 <- merge(mean.status,gams.data, by=c("Year","Site"))
m1 <- lmer(ihs(bio)~Microsite * gams + (1|Year), data=mean.spp)
shapiro.test(resid(m1))##
## Shapiro-Wilk normality test
##
## data: resid(m1)
## W = 0.9789, p-value = 0.8234anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 1.6259 1.6259 1 23.000 11.782 0.002272 **
## gams 18.1519 18.1519 1 23.367 131.533 4.448e-11 ***
## Microsite:gams 1.0755 1.0755 1 23.000 7.794 0.010365 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m1) ## R squared values ## R2m R2c
## 0.549095 0.945669m2 <- lmer(ihs(abd)~Microsite * gams + (1|Year), data=mean.spp)
shapiro.test(resid(m2))##
## Shapiro-Wilk normality test
##
## data: resid(m2)
## W = 0.97267, p-value = 0.6537car::Anova(m2, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ihs(abd)
## Chisq Df Pr(>Chisq)
## Microsite 0.1875 1 0.6650
## gams 63.2057 1 1.862e-15 ***
## Microsite:gams 0.1083 1 0.7421
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m2) ## R squared values ## R2m R2c
## 0.6343625 0.8521801m2.nat <- lmer(ihs(native)~Microsite * poly(gams,2) + (1|Year), data=mean.status2)
shapiro.test(resid(m2.nat))##
## Shapiro-Wilk normality test
##
## data: resid(m2.nat)
## W = 0.98084, p-value = 0.8701anova(m2.nat)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 1.2604 1.2604 1 22 1.1237 0.30063
## poly(gams, 2) 8.0514 4.0257 2 22 3.5890 0.04477 *
## Microsite:poly(gams, 2) 2.5294 1.2647 2 22 1.1275 0.34184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m2) ## R squared values ## R2m R2c
## 0.6343625 0.8521801m2.exo <- lmer(ihs(non.native)~Microsite * gams+ (1|Year), data=mean.status2)
shapiro.test(resid(m2.exo))##
## Shapiro-Wilk normality test
##
## data: resid(m2.exo)
## W = 0.90696, p-value = 0.01672anova(m2.exo)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 0.247 0.247 1 23.000 0.183 0.6727
## gams 74.389 74.389 1 23.965 55.056 1.175e-07 ***
## Microsite:gams 0.175 0.175 1 23.000 0.129 0.7223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m2) ## R squared values ## R2m R2c
## 0.6343625 0.8521801m3 <- lmer(rich~Microsite * poly(gams,2) + (1|Year), data=mean.spp)
shapiro.test(resid(m3))##
## Shapiro-Wilk normality test
##
## data: resid(m3)
## W = 0.93835, p-value = 0.1004car::Anova(m3, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: rich
## Chisq Df Pr(>Chisq)
## Microsite 0.6282 1 0.428012
## poly(gams, 2) 10.9223 2 0.004249 **
## Microsite:poly(gams, 2) 2.5197 2 0.283699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m3) ## R squared values ## R2m R2c
## 0.3445392 0.3495475## Merge MPD with gams
mpd.gams <- merge(mpd.rich, gams.data, by=c("Site","Year"))
m4 <- lmer(mpd~Microsite * gams + (1|Year), data=mpd.gams)
shapiro.test(resid(m4))##
## Shapiro-Wilk normality test
##
## data: resid(m4)
## W = 0.98022, p-value = 0.8555anova(m4)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 0.67216 0.67216 1 24 5.5363 0.02716 *
## gams 0.00313 0.00313 1 24 0.0258 0.87369
## Microsite:gams 0.54576 0.54576 1 24 4.4952 0.04452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m4) ## R squared values ## R2m R2c
## 0.2060199 0.2060199### Biomass
p1 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=ihs(bio)), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(bio)), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#E69F00", fill="#E69F0080", data=subset(mean.spp, Microsite=="open"), lwd=1.4) + stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#56B4E9", fill="#56B4E970", data=subset(mean.spp, Microsite=="shrub"), lwd=1.4) + ylab("Annnual plant biomass")+ xlab("") #annotate("text", x=-4.5,y=4, label="a", size=8)+coord_cartesian(ylim=c(0, 4))
### abundance
p2 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=ihs(abd)), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(abd)), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= ihs(y)~x,aes(x=gams, y=abd), color="black", fill="Grey50", data=mean.spp, lwd=1.4)+ ylab("Annual plant abundance")+ xlab("")+ annotate("text", x=-4.5,y=6, label="b", size=8)
### richness
p3 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=rich), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=rich), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#E69F00", fill="#E69F0080", data=subset(mean.spp, Microsite=="open"), lwd=1.4) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#56B4E9", fill="#56B4E970", data=subset(mean.spp, Microsite=="shrub"), lwd=1.4) + ylab("Species richness") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=4, label="c", size=8)
### phylogenetic diveristy
p4 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#56B4E9", size=3, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#E69F00", size=3, data=subset(mpd.gams, Microsite=="open")) + stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#56B4E9",fill="#56B4E980", lwd=1.4)+ stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#E69F00",fill="#E69F0080", lwd=1.4)+ ylab("Phylogenetic Community Dissimilarity") + xlab("Aridity Gradient")+ theme_Publication()+ annotate("text", x=-4.5,y=1.8, label="d", size=8)
grid.arrange(p1,p2,p3,p4)
Compare environmental differences - mixed model results
### Nutrients
nutrient.mean <- nutrients.climate %>% group_by(Site, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))
## potassium
m5 <- lm(potassium ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)##
## Call:
## lm(formula = potassium ~ poly(gams, 2) * microsite, data = nutrient.mean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.224 -15.052 -5.456 7.166 51.901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 141.46 12.13 11.662 2.67e-06 ***
## poly(gams, 2)1 -217.98 45.39 -4.803 0.001351 **
## poly(gams, 2)2 150.46 45.39 3.315 0.010617 *
## micrositeshrub 90.20 17.15 5.258 0.000766 ***
## poly(gams, 2)1:micrositeshrub -256.71 64.19 -4.000 0.003952 **
## poly(gams, 2)2:micrositeshrub 149.36 64.19 2.327 0.048389 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.09 on 8 degrees of freedom
## Multiple R-squared: 0.9641, Adjusted R-squared: 0.9416
## F-statistic: 42.95 on 5 and 8 DF, p-value: 1.438e-05anova(m5)## Analysis of Variance Table
##
## Response: potassium
## Df Sum Sq Mean Sq F value Pr(>F)
## poly(gams, 2) 2 170636 85318 82.835 4.502e-06 ***
## microsite 1 28476 28476 27.648 0.0007663 ***
## poly(gams, 2):microsite 2 22053 11026 10.706 0.0054741 **
## Residuals 8 8240 1030
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m5) ## R squared values ## R2m R2c
## 0.9429147 0.9429147p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=potassium), color="#56B4E9", size=3, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=potassium), color="#E69F00", size=3, data=subset(nutrient.mean, microsite=="open"))+ stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#56B4E9", fill="#56B4E980", data=subset(nutrient.mean, microsite=="shrub"), lwd=1.4) +
stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#E69F00", fill="#E69F0080",data=subset(nutrient.mean, microsite=="open"), lwd=1.4)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ annotate("text", x=-4.5,y=32, label="a", size=8)+coord_cartesian(ylim=c(0, 32))
se <- function(x) sd(na.omit(x))/sqrt(length(na.omit(x)))
swc <- smc %>% group_by(Year, Site, Microsite, Census) %>% summarize(swc.avg=mean(swc), swc.se=se(swc))
smc.gams <- merge(swc, gams.data, by=c("Site","Year"))
## early swc
smc.early <- subset(smc.gams, Census=="emergence")
m6 <- lmer(swc.avg~Microsite * poly(gams,2) + (1|Year), data=smc.early)
shapiro.test(resid(m6))##
## Shapiro-Wilk normality test
##
## data: resid(m6)
## W = 0.91432, p-value = 0.02516anova(m6)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 1.53 1.53 1 21.000 0.084 0.7752
## poly(gams, 2) 1449.15 724.58 2 21.343 39.560 6.617e-08 ***
## Microsite:poly(gams, 2) 3.69 1.85 2 21.000 0.101 0.9046
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m6) ## R squared values ## R2m R2c
## 0.5464095 0.8953768### soil moisture at emergence
p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#56B4E9", size=3, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#E69F00", size=3, data=subset(smc.early, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=1.4, color="black") + ylab("Soil moisture at emergence (%)") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=30, label="c", size=8)
## Create seasons for hobo
season1 <- subset(HOBO.data, Year==2015 & Month > 10 | Year==2016 & Month < 5)
season2 <- subset(HOBO.data, Year==2016 & Month > 10 | Year==2017 & Month < 5)
HOBO.season <- rbind(season1,season2)
HOBO.season[,"season"] <- c(rep("season.1",nrow(season1)),rep("season.2",nrow(season2)))
HOBO.season <- na.omit(HOBO.season)
hobo.means <- HOBO.season %>% group_by(season, Microsite, Site) %>% summarize(temp.avg=mean(Temp),temp.var=var(Temp),temp.sd=sd(Temp),temp.se=se(Temp),rh.avg=mean(RH, na.rm=T),rh.var=var(RH, na.rm=T),rh.se=se(RH),frost=(sum(na.omit(Temp)<0)/length(na.omit(Temp))*100),heat=(sum(na.omit(Temp)>30)/length(na.omit(Temp))*100))
hobo.means[,"Year"] <- c(rep(2016,14),rep(2017,14))
hobo.arid <- merge(hobo.means, gams.data, by=c("Site","Year"))
# hobo.arid[,"CV"] <-hobo.arid[,"temp.sd"] /hobo.arid[,"CV"]
## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))##
## Shapiro-Wilk normality test
##
## data: resid(m7)
## W = 0.97826, p-value = 0.8069anova(m7)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 1.51 1.51 1 21.000 0.255 0.6190
## poly(gams, 2) 707.64 353.82 2 21.104 59.704 2.051e-09 ***
## Microsite:poly(gams, 2) 0.11 0.05 2 21.000 0.009 0.9908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m7) ## R squared values ## R2m R2c
## 0.3627874 0.9560767p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#56B4E9", size=3, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#E69F00", size=3, data=subset(hobo.arid, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#E69F00", fill="#E69F0080", data=subset(hobo.arid, Microsite=="open"), lwd=1.4) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#56B4E9", fill="#56B4E970", data=subset(hobo.arid, Microsite=="shrub"), lwd=1.4) + ylab("Temperature variation") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=125, label="d", size=8)
shrubz <- read.csv("Data/ERG.shrub.csv")
shrub.mean <- shrubz %>% group_by(Site, Microsite) %>% summarize(compact=mean(Compaction))
shrub.clim <- data.frame(shrub.mean, gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))
#shrub.clim <- subset(shrub.clim, Year==2016)
m8 <- lm(log(compact) ~ Microsite * gams, data=shrub.clim)
summary(m8)##
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37348 -0.12596 0.04006 0.12843 0.33669
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.138777 0.804662 -0.172 0.8665
## Micrositeshrub -3.101869 1.137963 -2.726 0.0213 *
## gams 0.007791 0.010038 0.776 0.4556
## Micrositeshrub:gams 0.036028 0.014196 2.538 0.0295 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared: 0.6967, Adjusted R-squared: 0.6058
## F-statistic: 7.658 on 3 and 10 DF, p-value: 0.006anova(m8)## Analysis of Variance Table
##
## Response: log(compact)
## Df Sum Sq Mean Sq F value Pr(>F)
## Microsite 1 0.18836 0.18836 3.3176 0.098548 .
## gams 1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams 1 0.36569 0.36569 6.4409 0.029469 *
## Residuals 10 0.56776 0.05678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m8) ## R squared values ## R2m R2c
## 0.6386405 0.6386405p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=compact), color="#56B4E9", size=3, data=subset(shrub.clim, Microsite=="shrub"))+ geom_jitter(aes(x=gams, y=compact), color="#E69F00", size=3, data=subset(shrub.clim, Microsite=="open"))+ stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#56B4E9", fill="#56B4E980", data=subset(shrub.clim, Microsite=="shrub"), lwd=1.4) +
stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#E69F00", fill="#E69F0080",data=subset(shrub.clim, Microsite=="open"), lwd=1.4)+ ylab("Soil compaction") + xlab("")+ theme_Publication() + annotate("text", x=-4.4,y=2.5, label="b", size=8) +coord_cartesian(ylim=c(0.5, 2.5))
grid.arrange(p5,p8,p6,p7)
### Microsite comparisons
# nut <- read.csv("Data//ERG.soilnutrients.csv")
#
# fit1 <- aovp(N ~ microsite, data=nut)
# summary(fit1)
# t.test(N ~ microsite, data=nut)
# fit2 <- aovp(P ~ microsite, data=nut)
# t.test(P ~ microsite, data=nut)
# summary(fit2)
# fit3 <- aovp(K ~ microsite, data=nut)
# t.test(K ~ microsite, data=nut)
# summary(fit3)
#
# compact <- read.csv("Data//ERG.shrub.csv")
#
# fit4 <- aovp(Compaction ~ Microsite, data=compact)
# t.test(Compaction ~ Microsite, data=compact)
# summary(fit4)
#
# soilmoist <- read.csv("Data//ERG.phytometer.census.csv")
#
# fit5 <- aovp(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
# summary(fit5)
# t.test(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
#
# micro <- read.csv("Data//ERG.logger.data.csv")
#
# micro.avg <- micro %>% group_by(Year, Site, Microsite, Rep) %>% summarize(temp=mean(Temp), rh=mean(RH))
#
# fit6 <- aovp(temp ~ Microsite, data=micro.avg)
# summary(fit6)
# t.test(temp ~ Microsite, data=micro.avg)
#
# fit7 <- aovp(rh ~ Microsite, data=micro.avg)
# summary(fit7)
# t.test(rh ~ Microsite, data=micro.avg)Plot level GLMM
m1 <- glmer.nb(Abundance ~ Microsite * gams + (1|Year), data=community.arid)
car::Anova(m1, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Abundance
## Chisq Df Pr(>Chisq)
## Microsite 16.3613 1 5.234e-05 ***
## gams 1136.3911 1 < 2.2e-16 ***
## Microsite:gams 1.0278 1 0.3107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coef(m1)## $Year
## (Intercept) Micrositeshrub gams Micrositeshrub:gams
## 2016 13.13187 0.6544052 -0.1298036 -0.005782823
## 2017 11.53434 0.6544052 -0.1298036 -0.005782823
##
## attr(,"class")
## [1] "coef.mer"source("rsquared.r")
m2 <- glmer.nb(Richness ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
car::Anova(m2, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Richness
## Chisq Df Pr(>Chisq)
## Microsite 6.31 1 0.01201 *
## poly(gams, 2) 1036.41 2 < 2e-16 ***
## Microsite:poly(gams, 2) 374.26 2 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rsquared.glmm(m2)
# m2.gam <- gam(Richness~ Microsite +s(gams, k=3) + s(gams, by=Microsite, k=3) + s(Year, bs="re"), family=poisson, data=community.arid)
p1 <- ggplot(community.arid, aes(x=gams, y=Richness)) + geom_jitter(aes(x=gams, y=Richness), color="#181818", size=3, alpha=0.4, width = 0.75,shape = 1, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=Richness), color="#181818", size=2,alpha=0.4, data=subset(community.arid, Microsite=="open"),shape = 2, width = 0.75) + theme_Publication() +
stat_smooth(method="gam", method.args=list(family="poisson"), formula= y~s,aes(x=gams, y=Richness), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)
## richness
p1 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=Richness), color="#181818", size=3, alpha=0.4, width = 0.75,shape = 1, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=Richness), color="#181818", size=2,alpha=0.4, data=subset(community.arid, Microsite=="open"),shape = 2, width = 0.75) + theme_Publication() + ylab("species richness") + xlab("Gams index of continentality")+
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080" , data=subset(community.arid, Microsite=="shrub"), lwd=2) +
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+ ylab("Species richness") + xlab("") + annotate("text", x=50,y=6, label="a", size=8)
m3<- lmer(ihs(Biomass) ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
anova(m3)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 64.53 64.527 1 833.00 210.49 < 2.2e-16 ***
## poly(gams, 2) 511.61 255.805 2 833.44 834.46 < 2.2e-16 ***
## Microsite:poly(gams, 2) 32.08 16.041 2 833.00 52.33 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m3) ## R squared values## R2m R2c
## 0.5031296 0.8793048# m3.gam <- gam(ihs(Biomass)~ s(gams, by=Microsite, k=3) + s(Year, bs="re"), family=gaussian, data=community.arid)
#
## biomass
p2 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=3, width = 0.75,shape = 1, alpha=0.5, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=2,shape = 2,alpha=0.5, width = 0.75, data=subset(community.arid, Microsite=="open")) + theme_Publication() +
geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="shrub"), lwd=2) +
geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+ ylab("Annual biomass") + xlab("")+ annotate("text", x=50,y=5, label="b", size=8)
m4<- lmer(mpd~ Microsite * gams + (1|Year), data=mpd.gams)
anova(m4)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 0.67216 0.67216 1 24 5.5363 0.02716 *
## gams 0.00313 0.00313 1 24 0.0258 0.87369
## Microsite:gams 0.54576 0.54576 1 24 4.4952 0.04452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m4) ## R squared values## R2m R2c
## 0.2060199 0.2060199## phylogenetic diversity
p3 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#181818", size=3,shape=1, alpha=0.4, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#181818",alpha=0.4, size=2,shape=2, data=subset(mpd.gams, Microsite=="open")) + stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#181818",fill="#80808080", lwd=2)+ stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#181818",fill="#80808080", lwd=2, lty=2)+ ylab("Phylogenetic community dissimilarity") + xlab("Gams index of continentality")+ theme_Publication()+ annotate("text", x=50,y=1.4, label="c", size=8)
## native vs non-native
status.gams <- merge(status.data, gams.data, by=c("Year","Site"))
m5<- glmer.nb(native~ Microsite * scale(gams) + (1|Year), data=status.gams) ## need to scale data
car::Anova(m5, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: native
## Chisq Df Pr(>Chisq)
## (Intercept) 565.724 1 < 2.2e-16 ***
## Microsite 23.122 1 1.521e-06 ***
## scale(gams) 67.393 1 2.224e-16 ***
## Microsite:scale(gams) 35.708 1 2.292e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rsquared.glmm(m5)
m6<- glmer.nb(non.native~ Microsite * scale(gams) + (1|Year), data=status.gams) ## need to scale data
car::Anova(m6, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: non.native
## Chisq Df Pr(>Chisq)
## (Intercept) 12.5323 1 0.000400 ***
## Microsite 13.5445 1 0.000233 ***
## scale(gams) 353.8669 1 < 2.2e-16 ***
## Microsite:scale(gams) 1.4368 1 0.230664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rsquared.glmm(m5)
status.plot <- gather(status.gams, origin, abundance, 5:6)
p4 <- ggplot(subset(status.plot, abundance>0)) + geom_jitter(aes(x=gams, y=abundance), color="#181818", size=3, data=subset(status.plot, Microsite=="shrub"), width = 1, shape=1, alpha=0.3) +geom_jitter(aes(x=gams, y=abundance), width = 1, color="#181818", size=2,alpha=0.3,shape=2, data=subset(status.plot, Microsite=="open")) + theme_Publication() +
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance), color="#505050", se=F, data=subset(status.plot, Microsite=="open" & origin=="native"), lwd=2, lty=2)+
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818", data=subset(status.plot, Microsite=="open" & origin=="non.native"), lwd=2, lty=2)+
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#505050", data=subset(status.plot, Microsite=="shrub" & origin=="native"), lwd=2)+
geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818", data=subset(status.plot, Microsite=="shrub" & origin=="non.native"), lwd=2) + ylab("Annual plant abundance") + xlab("Gams index of continentality")+ annotate("text", x=50,y=310, label="d", size=8)
grid.arrange(p1,p2,p3,p4)
### Nutrients
nutrient.mean <- nutrients.climate %>% group_by(Site, gams, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))
## nitrogen
m5 <- lm(nitrogen ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)##
## Call:
## lm(formula = nitrogen ~ poly(gams, 2) * microsite, data = nutrient.mean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2918 -0.9251 -0.2773 0.8012 5.6374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.595 1.038 1.536 0.163115
## poly(gams, 2)1 3.927 3.885 1.011 0.341676
## poly(gams, 2)2 3.847 3.885 0.990 0.350969
## micrositeshrub 7.676 1.468 5.228 0.000795 ***
## poly(gams, 2)1:micrositeshrub 22.258 5.494 4.052 0.003676 **
## poly(gams, 2)2:micrositeshrub 17.833 5.494 3.246 0.011771 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.747 on 8 degrees of freedom
## Multiple R-squared: 0.9298, Adjusted R-squared: 0.8859
## F-statistic: 21.18 on 5 and 8 DF, p-value: 0.0002012anova(m5)## Analysis of Variance Table
##
## Response: nitrogen
## Df Sum Sq Mean Sq F value Pr(>F)
## poly(gams, 2) 2 389.59 194.794 25.817 0.0003239 ***
## microsite 1 206.22 206.223 27.332 0.0007948 ***
## poly(gams, 2):microsite 2 203.35 101.677 13.476 0.0027446 **
## Residuals 8 60.36 7.545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m5) ## R squared values ## R2m R2c
## 0.8906817 0.8906817p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3, shape=16,alpha=0.5, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3,shape=17, alpha=0.4, data=subset(nutrient.mean, microsite=="open"))+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080", data=subset(nutrient.mean, microsite=="shrub"), lwd=2) +
stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080",data=subset(nutrient.mean, microsite=="open"), lwd=2, lty=2)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ coord_cartesian(ylim=c(-2, 26)) + xlab("")+ annotate("text", x=63,y=25, label="a", size=8)
smc.early <- subset(smc, Census=="emergence")
smc.early <- merge(smc.early, gams.data, by=c("Year","Site"))
smc.early <- smc.early %>% group_by(Year, Site, gams, Microsite) %>% summarize(swc.avg=mean(swc))
## nitrogen
m6 <- lmer(log(swc.avg) ~ poly(gams,2)* Microsite + (1|Year), data=smc.early)
summary(m6)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: log(swc.avg) ~ poly(gams, 2) * Microsite + (1 | Year)
## Data: smc.early
##
## REML criterion at convergence: 16.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75790 -0.67057 0.05989 0.64573 1.52776
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year (Intercept) 0.32730 0.5721
## Residual 0.09642 0.3105
## Number of obs: 28, groups: Year, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 2.183930 0.412962 1.000000 5.288
## poly(gams, 2)1 -4.095881 0.529751 21.483000 -7.732
## poly(gams, 2)2 0.009355 0.444206 21.040000 0.021
## Micrositeshrub -0.021937 0.117363 21.000000 -0.187
## poly(gams, 2)1:Micrositeshrub 0.108296 0.621024 21.000000 0.174
## poly(gams, 2)2:Micrositeshrub -0.128193 0.621024 21.000000 -0.206
## Pr(>|t|)
## (Intercept) 0.119
## poly(gams, 2)1 1.22e-07 ***
## poly(gams, 2)2 0.983
## Micrositeshrub 0.854
## poly(gams, 2)1:Micrositeshrub 0.863
## poly(gams, 2)2:Micrositeshrub 0.838
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(,2)1 pl(,2)2 Mcrsts p(,2)1:
## ply(gms,2)1 0.000
## ply(gms,2)2 0.000 0.084
## Microstshrb -0.142 0.000 0.000
## ply(g,2)1:M 0.000 -0.586 0.000 0.000
## ply(g,2)2:M 0.000 0.000 -0.699 0.000 0.000anova(m6)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## poly(gams, 2) 8.6913 4.3456 2 21.339 45.071 2.183e-08 ***
## Microsite 0.0034 0.0034 1 21.000 0.035 0.8535
## poly(gams, 2):Microsite 0.0070 0.0035 2 21.000 0.037 0.9642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m6) ## R squared values ## R2m R2c
## 0.5883174 0.9063213p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3, shape=16,alpha=0.5, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3,shape=17,alpha=0.5, data=subset(smc.early, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=2, color="black") + ylab("Soil moisture at emergence (%)") + xlab("")+ annotate("text", x=50,y=35, label="b", size=8)
## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))##
## Shapiro-Wilk normality test
##
## data: resid(m7)
## W = 0.97826, p-value = 0.8069anova(m7) ## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Microsite 1.51 1.51 1 21.000 0.255 0.6190
## poly(gams, 2) 707.64 353.82 2 21.104 59.704 2.051e-09 ***
## Microsite:poly(gams, 2) 0.11 0.05 2 21.000 0.009 0.9908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m7) ## R squared values ## R2m R2c
## 0.3627874 0.9560767p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=16, alpha=0.5, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=17, alpha=0.5, data=subset(hobo.arid, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="open"), lwd=2, lty=2) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="shrub"), lwd=2) + ylab("Temperature variation") + xlab("Gams index of continentality")+ annotate("text", x=50,y=125, label="c", size=8)
## soil compaction
m8 <- lm(log(compact) ~ Microsite * gams, data=shrub.clim)
summary(m8)##
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37348 -0.12596 0.04006 0.12843 0.33669
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.138777 0.804662 -0.172 0.8665
## Micrositeshrub -3.101869 1.137963 -2.726 0.0213 *
## gams 0.007791 0.010038 0.776 0.4556
## Micrositeshrub:gams 0.036028 0.014196 2.538 0.0295 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared: 0.6967, Adjusted R-squared: 0.6058
## F-statistic: 7.658 on 3 and 10 DF, p-value: 0.006anova(m8)## Analysis of Variance Table
##
## Response: log(compact)
## Df Sum Sq Mean Sq F value Pr(>F)
## Microsite 1 0.18836 0.18836 3.3176 0.098548 .
## gams 1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams 1 0.36569 0.36569 6.4409 0.029469 *
## Residuals 10 0.56776 0.05678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m8) ## R squared values ## R2m R2c
## 0.6386405 0.6386405p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=16, alpha=0.5, data=subset(shrub.clim, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=17, alpha=0.5, data=subset(shrub.clim, Microsite=="open")) + theme_Publication()+ stat_smooth(method="lm", formula= y~x, aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="open"), lwd=2, lty=2) + stat_smooth(method="lm", formula= y~x,aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="shrub"), lwd=2) + ylab("Soil compaction") + xlab("Gams index of continentality")+ annotate("text", x=63,y=1.1, label="d", size=8)
grid.arrange(p5,p6,p7,p8)
Phytometer with GAMS
census <- read.csv("Data/ERG.phytometer.census.csv")
census[is.na(census)] <- 0
census.arid <- merge(census,gams.data, by=c("Year","Site"))
## treatment plots
census.long <- gather(census.arid, species, biomass, 12:14)
census.long[,"species"] <- ifelse(census.long[,"species"] =="Phacelia.biomass", "P. tanacetifolia",
ifelse(census.long[,"species"]=="Plantago.biomass","P. insularis","S. columbariae"))
census.long[,"species"] <- as.factor(census.long$species)
## microsite average
census.micro <- census.long %>% filter(biomass>0) %>% group_by(Microsite, species) %>% summarize(Biomass=mean(biomass), error=se(biomass))
p9 <- ggplot(census.micro, aes(x=species, y=Biomass, fill=Microsite)) + geom_bar(position=position_dodge(), stat="identity", colour="black") + scale_fill_brewer(palette="Greys")+ geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
width=.2, # Width of the error bars
position=position_dodge(.9))+ theme_Publication() + theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## nutrient average
census.nut <- census.long %>% filter(biomass>0) %>% group_by(Nutrient, species) %>% summarize(Biomass=mean(biomass), error=se(biomass))
p10 <- ggplot(census.nut, aes(x=species, y=Biomass, fill=Nutrient)) + geom_bar(position=position_dodge(), stat="identity",colour="black") + scale_fill_brewer(palette="Greys")+ geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
width=.2, # Width of the error bars
position=position_dodge(.9))+ theme_Publication()+ theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## logistic first
census.long[,"presence"] <- ifelse(census.long$biomass>0,1,0)
## Drop panoche 2017 because of cold
census.npan <- census.long %>% filter(Site!="PanocheHills" | Year!="2017")
m1.l <- glmer(presence~Microsite * gams * Nutrient + (1|Year),data=subset(census.npan, species=="P. tanacetifolia"), family=binomial)
car::Anova(m1.l, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: presence
## Chisq Df Pr(>Chisq)
## Microsite 23.7693 1 1.086e-06 ***
## gams 6.7832 1 0.009202 **
## Nutrient 0.0047 1 0.945204
## Microsite:gams 3.2156 1 0.072938 .
## Microsite:Nutrient 1.6347 1 0.201055
## gams:Nutrient 0.0000 1 0.995620
## Microsite:gams:Nutrient 0.0062 1 0.937144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. tanacetifolia" & presence==1))
shapiro.test(resid(m1))##
## Shapiro-Wilk normality test
##
## data: resid(m1)
## W = 0.99225, p-value = 0.2595anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value
## Microsite 22.121 22.121 1 221 12.213
## poly(gams, 2) 157.099 78.549 2 221 43.368
## Nutrient 5.710 5.710 1 221 3.152
## Microsite:poly(gams, 2) 15.688 7.844 2 221 4.331
## Microsite:Nutrient 0.078 0.078 1 221 0.043
## poly(gams, 2):Nutrient 2.854 1.427 2 221 0.788
## Microsite:poly(gams, 2):Nutrient 1.406 0.703 2 221 0.388
## Pr(>F)
## Microsite 0.0005733 ***
## poly(gams, 2) < 2.2e-16 ***
## Nutrient 0.0771914 .
## Microsite:poly(gams, 2) 0.0142925 *
## Microsite:Nutrient 0.8361223
## poly(gams, 2):Nutrient 0.4561215
## Microsite:poly(gams, 2):Nutrient 0.6786980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m1) ## R squared values ## R2m R2c
## 0.346015 0.346015m2.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="P. insularis"), family=binomial)
car::Anova(m2.l, type=2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: presence
## Chisq Df Pr(>Chisq)
## Microsite 1.4664 1 0.22591
## poly(gams, 2) 8.3747 2 0.01519 *
## Nutrient 0.2723 1 0.60178
## Microsite:poly(gams, 2) 4.8534 2 0.08833 .
## Microsite:Nutrient 1.6652 1 0.19691
## poly(gams, 2):Nutrient 2.9034 2 0.23418
## Microsite:poly(gams, 2):Nutrient 0.4909 2 0.78236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="P. insularis" & presence==1))
shapiro.test(resid(m2))##
## Shapiro-Wilk normality test
##
## data: resid(m2)
## W = 0.99193, p-value = 0.5364anova(m2)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value
## Microsite 0.569 0.5693 1 141.23 0.5147
## poly(gams, 2) 36.153 18.0765 2 137.31 16.3429
## Nutrient 8.958 8.9576 1 141.27 8.0986
## Microsite:poly(gams, 2) 0.858 0.4289 2 141.51 0.3878
## Microsite:Nutrient 0.110 0.1096 1 141.02 0.0991
## poly(gams, 2):Nutrient 2.373 1.1866 2 141.20 1.0728
## Microsite:poly(gams, 2):Nutrient 0.654 0.3272 2 141.03 0.2958
## Pr(>F)
## Microsite 0.47430
## poly(gams, 2) 4.297e-07 ***
## Nutrient 0.00509 **
## Microsite:poly(gams, 2) 0.67929
## Microsite:Nutrient 0.75338
## poly(gams, 2):Nutrient 0.34481
## Microsite:poly(gams, 2):Nutrient 0.74439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m2) ## R squared values ## R2m R2c
## 0.2802814 0.5695941m3.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="S. columbariae"), family=binomial)
car::Anova(m3.l, type=3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: presence
## Chisq Df Pr(>Chisq)
## (Intercept) 55.9936 1 7.271e-14 ***
## Microsite 1.2916 1 0.2558
## poly(gams, 2) 18.5593 2 9.330e-05 ***
## Nutrient 0.0728 1 0.7873
## Microsite:poly(gams, 2) 1.7247 2 0.4222
## Microsite:Nutrient 2.3865 1 0.1224
## poly(gams, 2):Nutrient 0.8524 2 0.6530
## Microsite:poly(gams, 2):Nutrient 0.5872 2 0.7456
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m3 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="S. columbariae" & biomass>0))
shapiro.test(resid(m3))##
## Shapiro-Wilk normality test
##
## data: resid(m3)
## W = 0.99363, p-value = 0.1532anova(m3)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value
## Microsite 0.350 0.3503 1 333.03 0.2070
## poly(gams, 2) 61.028 30.5139 2 327.98 18.0272
## Nutrient 18.797 18.7973 1 333.12 11.1052
## Microsite:poly(gams, 2) 1.304 0.6521 2 333.03 0.3853
## Microsite:Nutrient 0.386 0.3859 1 333.05 0.2280
## poly(gams, 2):Nutrient 2.914 1.4569 2 333.05 0.8607
## Microsite:poly(gams, 2):Nutrient 3.356 1.6780 2 333.01 0.9913
## Pr(>F)
## Microsite 0.6494522
## poly(gams, 2) 3.733e-08 ***
## Nutrient 0.0009576 ***
## Microsite:poly(gams, 2) 0.6805772
## Microsite:Nutrient 0.6333374
## poly(gams, 2):Nutrient 0.4237872
## Microsite:poly(gams, 2):Nutrient 0.3721737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r.squaredGLMM(m3) ## R squared values ## R2m R2c
## 0.1506253 0.2337410## summary presence
mean.occ <- census.npan %>% group_by(gams,species, Microsite) %>% summarize(presence=mean(presence))
p11 <- ggplot(census.npan, aes(x=gams, y=presence, fill=species, color=species)) + theme_Publication()+
geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2", data=subset(census.npan, species=="P. tanacetifolia" & Microsite=="shrub"), lwd=2)+
geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2", data=subset(census.npan, species=="P. tanacetifolia"& Microsite=="open"), lwd=2, lty=2) + geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence), color="#fb6f70", fill="#fb6f7050", data=subset(census.npan, species=="P. insularis"), lwd=2)+ geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence),fill="#b7d73850", color="#b7d738", data=subset(census.npan, species=="S. columbariae"), lwd=2)+ ylab("Probability of occurrence") + xlab("Gams index of continentality")
p12 <- ggplot(census.long, aes(x=gams, y=biomass, fill=species, color=species)) + theme_Publication() + geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)), color="#fb6f70", fill="#fb6f7050", data=subset(census.long, species=="P. insularis" & biomass>0), lwd=2, fullrange=TRUE) + geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#b7d73850", color="#b7d738", data=subset(census.long, species=="S. columbariae"& biomass>0), lwd=3)+
geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2", data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="shrub"), lwd=2, fullrange=TRUE)+geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2", data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="open"), lwd=2,lty=2, fullrange=TRUE)+ ylab("Phytometer biomass") + xlab("Gams index of continentality")
grid.arrange(p9,p11,p10, p12)